//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
// g++ -I /home/tiago/astro_local/lib/ -O2 cur_luz.cpp -o cur_luz
//
// cur_luz.cpp - Nova versao do programa curluz2.c agora escrito em C++.
// Inclui o biblioteca de leitura de arquivos dataFileReader.h (Zabot A.)
// de modo que comentarios podem ser inseridos nos arquivos de entrada. Os
// comentarios devem ser precedidos de #.
// 
// Modo de utilizacao do programa:
// cur_luz arq_in arq_out
// onde arq_in o arquivo com os parametros de entrada e arq_out o arquivo 
// para serem escritos os resultados.
//
// Formato do arquivo de entrada:
// 1a. linha de dados (sem contar linha de comentario):
//   i  - inclinacao (graus).
//   q  - razao de massa (M2/M1).
//   f2 - fluxo da secundaria em fase 0. 
//   fd - fluxo do disco (constante somada a curva de luz).
//   S/N- Relacao Sinal-Ruido na curva
// 2a. linha de dados:
//   limb_law  - lei de obscurecimento de borda.
//   grav_coef - coeficiente do obscurecimento gravitacional.
//   Tpole     - Temperatura do polo da estrela.
//   nspot     - Numero de spots 
//   f_fac     -
// 3a. linha de dados:
//   (coeficientes de obscurecimento de borda)
// 4a. a (4+nspot) linhas de dados: 
//   theta - posicao theta do spot
//   phi   - posicao phi do spot
//   dt    - largura a meia altura em theta
//   dp    - largura a meia altura em phi
//   Is    - Intensidade em relacao a I_L1
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+

#include <iostream>
#include <cmath>
#include <fstream>
#include <cmath>
#include <vector>
#include <sstream>
#include <string>
#include <cstring>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_blas.h>
#include "light_curve.h"
//#include "light_curve.cpp"
//#include "roche.cpp"
#include "convolve.hpp"

#include "dataFileReader_v0.1.1/dataFileReader.h"

//#include "lib_zabot/dataFileReader.h"


using namespace std;

int definicoes_obscurecimento(const int law);
void show_par(ostream& out, vector<vector<double> > &par);
int find_io(int argc, char *argv[]);
void inline help_and_exit(void);
 
int main(int argc, char *argv[])
{
	int stat = setenv("GSL_RNG_SEED",ntos<int>(getpid()).c_str(),1);
	
	char * descr = getenv("GSL_RNG_SEED");
	if (descr)
		cout<<"# - Value of GSL_RNG_SEED is: "<< descr<<endl;
	else {
		cout<<"# - GSL_RNG_SEED not defined, using default"<<endl;
	}
	
  // Testar argumentos de entrada...
  // Definindo qual a forma de input-output! Mostra help e exit() caso indefinido...
  // 1 - arq_in arq_out
  // 2 - arq_in arq_out fases.dat
  // 3 - arq_in arq_out image.img
  // 4 - arq_in arq_out fases.dar image.img

  int IO_type = find_io(argc, argv);

  cout << "# Modo In - Out: " << IO_type << endl;

  // Input File Definition - IFD:  
  //     Number of fixed Parameters lines      -  IFD_NP
  //     Number of columns Parameters          - *IFD_P (array[IFD_NP])
  //     Spot Line(onde com. def. de spot - 1) -  IFD_SL
  //     Number of Spot Parameters             -  IFD_SP
  //     posicao da lei de obs. de Borda       -  IFD_B

  const int IFD_NP = 2;
  const int IFD_L1 = 1; 
  const int IFD_P[IFD_NP] = {9, 7};
  const int IFD_SL = 3;
  const int IFD_SP = 5;
  const int IFD_B = 2;
  int SIZE = 4;

  // Input Image Definition - IID: Sobre imagens de saida de roch_tom.cpp

  const int IID_NC = 12; // Numero de Colunas

  DataFileReader fin(argv[1]);
  ++fin;
  vector<vector<double> > in_par(IFD_NP+1);

  //Lendo resolucao do mapa

  if(fin.numberOfColumns() == IFD_L1){
    //fin >> SIZE;
	  SIZE = fin.readColumn<int>(1);
    ++fin;
  }
  else
    cout << "# - Utilizando tamanho de mapa padrao...\n";
  
  cout << "# - SIZE = " << SIZE << endl;

  //Lendo parametros orbitais

  for(int i = 0; i < IFD_NP; i++){
    if(fin.numberOfColumns() == IFD_P[i]){
      double data_in;
      for(int j = 0; j < IFD_P[i]; j++){
	//fin >> data_in;
		  data_in = fin.readColumn<double>(j+1);
	if(i == 0 && j == 0) data_in *= M_PI/180.;
	in_par[i].push_back(data_in);
      }
    }
    else { 
      ostringstream erro;
      erro << "Numero de parametros errado\n\tArquivo: " 
	   << argv[1] << "\n\tLinha  : " << fin.lineNumber();
      aborta(erro.str());
    }
    ++fin;
/*    if(fin.eof()) {
      ostringstream erro;
      erro << "Numero de parametros errado\n\tArquivo: " 
	   << argv[1] << "\n\tLinha  : " << fin.lineNumber();
      aborta(erro.str());
    }*/
  }

  //Lendo coeficientes do obscurecimento de borda

  if(fin.numberOfColumns()==definicoes_obscurecimento(int(in_par[1][IFD_B]))){
    double data_in;
    for(int i = 0; i < definicoes_obscurecimento(int(in_par[1][IFD_B])); i++){
      //fin >> data_in;
		data_in = fin.readColumn<double>(i+1);
      in_par[2].push_back(data_in);
    }
  }
  else {
    ostringstream erro;
    erro << "Numero de parametros errado\n\tArquivo: " 
	 << argv[1] << "\n\tLinha  : " << fin.lineNumber();
    aborta(erro.str());
  }

  ++fin;
  /*
  if(fin.eof()) {
    ostringstream erro;
    erro << "Numero de parametros errado\n\tArquivo: " 
	 << argv[1] << "\n\tLinha  : " << fin.lineNumber();
    aborta(erro.str());
  }
  */
  //Lendo parametros do spot

  for(int i = 0; i < int(in_par[1][3]); i++){
    if(fin.numberOfColumns() == IFD_SP){
      double data_in;
      in_par.push_back(vector<double> ());
      for(int j = 0; j < IFD_SP; j++){
	//fin >> data_in;
		  data_in = fin.readColumn<double>(j+1);
	in_par[IFD_SL+i].push_back(data_in);
      }
    }
    else {
      ostringstream erro;
      erro << "Numero de parametros errado\n\tArquivo: " 
	   << argv[1] << "\n\tLinha  : " << fin.lineNumber();
      aborta(erro.str());
    }
    ++fin;
  }

  fin.close();

  show_par(cout, in_par);

  if(in_par[0][4] == 0) in_par[0][4] == 1; //se sinal ruido = 0 S/R = 1

  vector<double> lin_form(3);
  lin_form[0] = -in_par[0][5];
  lin_form[1] = in_par[0][5];
  lin_form[2] = in_par[0][6];

  const gsl_rng_type * T;
  gsl_rng * r;

  gsl_rng_env_setup();

  T = gsl_rng_default;
  r = gsl_rng_alloc (T);

  if(in_par[0][4] > 0) {
    double sn = in_par[0][4]; 
    in_par[0][4] = (1./sn);
  }

  //gaussiana para blurring instrumental
  vector<vector<double> > f2(2);
  f2[0].resize(floor(in_par[0][5]*200./in_par[0][6]));
  f2[1].resize(floor(in_par[0][5]*200./in_par[0][6]));
  //double dxx = (in_par[0][5]*2.0/in_par[0][6]);
  //cout << dxx << endl;
  int itr=0;
  for ( double x0 = -in_par[0][5]; x0 < in_par[0][5]; x0+=in_par[0][6]/100.){
    f2[0][itr] = x0;
    f2[1][itr] = exp(-(x0/in_par[0][7])*(x0/in_par[0][7]));
    itr++;
    if(itr > f2[0].size()) break;
  }
	cout << itr << " " << f2[0].size() <<endl;
//  for(int i = 0; i < f2[0].size(); i++)
//    cout << f2[0][i] << " " << f2[1][i] << endl;
//  exit(0);
  switch (IO_type) {
  case 1:
    {
    Roche cl2(SIZE, in_par);
    cl2.vorb(false);
    cl2.refaz_super();
    //se tiver spot manda fazer!
    if(in_par[1][3] > 0) cl2.spot_maker();

    for(double fase = 0; fase <= 1; fase += 0.05){
      //prepara e abre arquivo de saida
      ostringstream saida;
      saida << argv[2] << '_' << fase;
		double flux2=cl2.flux(fase);
      //cout << fase << "  " << flux2 << endl;
      ofstream oo(saida.str().c_str());
      
      if(!oo) {
	ostringstream erro;
	erro << "Arquivo de saida nao pode ser aberto\n\tArquivo: " 
	     << argv[2];
	aborta(erro.str());
      }

      //escreve saida no arquivo
      
      show_par(oo, in_par);
      
      Histograma lp = cl2.rad_prof(fase, lin_form); //lin_prof;
      vector<vector<double> > line_mod(2);
      line_mod[0].resize(lp.size());
      line_mod[1].resize(lp.size());
      //cout << "#- Aqui 0\n";
      for(int i = 0; i < lp.size(); i++){
        line_mod[0][i] = lp.le1(i);
        line_mod[1][i] = -lp.le2(i)+flux2+in_par[1][6];
      }
      //cout << "#- Aqui 1\n";
      //tmp_vec=vector<double> (lp.col2());
      //line_mod[1].push_back(tmp_vec);
      //se tiver blurring aplica...
      if(in_par[0][7] > 0.0){
        Convolve ff(line_mod,f2,line_mod[0].size(),f2[0].size());
        ff.xmin(line_mod[0][0]);
        ff.xmax(line_mod[0][line_mod[0].size()-1]);
        for (int i = 0; i < line_mod[0].size(); i++){
          gsl_function F;
          F.function = &cnv_fnc;
          F.params = &ff;
          GSL_Integration itr_conv(F);
          itr_conv.xmin(ff.xmin());
          itr_conv.xmax(ff.xmax());
          ff.px(line_mod[0][i]);
          line_mod[1][i] = itr_conv.integrate_qag()/in_par[0][7]/sqrt(M_PI);
          itr_conv.clear();
          //cout << line_mod[0][i]
          //<< " " << line_mod[1][i] << " " 
          //<< itr_conv.integrate_qag()/2.0 << endl;
        }
	
      }
      cout << "#Escrevendo saida em " << saida.str() << " ...\n";

      for (int i = 2; i < lp.size()-2; i++){
	int iter = (i != (lp.size()-1) ? 1 : -1);
	//if(lp.le2(i) > 0 || lp.le2(i+iter) > 0) 
		  oo << line_mod[0][i]/1.0e5<< "  " << line_mod[1][i]+((line_mod[1][i])*in_par[0][4]*(2.*gsl_rng_uniform(r) - 1. ))
		  << " " << (line_mod[1][i])*in_par[0][4] << endl;

 // oo << 5000.0*(1+line_mod[0][i]/299792458.0e2) << "  " << line_mod[1][i]+((line_mod[1][i])*in_par[0][4]*(2.*gsl_rng_uniform(r) - 1. ))
//	     << " " << (line_mod[1][i])*in_par[0][4] << endl;
	/*oo << line_mod[0][i]
	     << "  " << line_mod[1][i]+((line_mod[1][i])*in_par[0][4]*(2.*gsl_rng_uniform(r) - 1. ))
	     << " " << (line_mod[1][i])*in_par[0][4] << endl;*/
      }
      //fecha saida
      oo.close();
    }
    }
    break;
  case 2: //caso de arquivo com fases dados
    {
    Roche cl2(SIZE, in_par);
    double c = 3e10; // cm/s
    double y0 = 5000.0; // Angstrons
    //cl2.vorb(false);
	cl2.vorb(true);
    cl2.refaz_super();
	cl2.CalcPixSize();
	cout << "# - Npix = " << cl2.npix() << endl
		<< "# - Apix = " << cl2.pix_size() << endl;
    //se tiver spot manda fazer!
    if(in_par[1][3] > 0) {
		cl2.spot_maker();
		/*for (int ispot = 0; ispot<in_par[1][3]; ispot++) {
			cout << in_par[3+ispot][0] << " " 
			<< in_par[3+ispot][1] << " "
			<< in_par[3+ispot][4] << endl ;
			cl2.sum_flux_pix(cl2.proc_pix(in_par[3+ispot][0],in_par[3+ispot][1]),in_par[3+ispot][4]);
		}*/
		
	}

    DataFileReader f_fase(argv[3]);
    while(!++f_fase){
      //leio fase do arquivo
      double fase;
      //f_fase >> fase;
		fase = f_fase.readColumn<double>(1);
      double flux2=cl2.flux(fase);
      //preparo saida
      ostringstream saida;
      saida << argv[2] << '_' << fase;
      //cout << fase << "  " << flux2 << endl;
      ofstream oo(saida.str().c_str());
      
      if(!oo) {
	ostringstream erro;
	erro << "Arquivo de saida nao pode ser aberto\n\tArquivo: " 
	     << argv[2];
	aborta(erro.str());
      }

      //escreve saida no arquivo

		show_par(oo, in_par);
		
		
		int _Np_ = ceil(in_par[0][8] / ( in_par[0][6]/1e5/(2.*M_PI*cl2.K2()) ) ); //Numero de perfis para combinar
		
		if (_Np_ < 2) { //This is as lower as it can go!
			_Np_ = 1;
		}
		cout << "# - Np = "<< _Np_ << endl;
		
		Histograma* lp = new Histograma [_Np_]; //cl2.rad_prof(fase, lin_form); //lin_prof;

		for (int nhist = 0; nhist < _Np_; ++nhist) {
			double dphi = _Np_ > 1 ? in_par[0][8]/(_Np_-1) * nhist - in_par[0][8]/2.0 : 0.0;

			lp[nhist] = cl2.rad_prof(fase + dphi , lin_form);
		}

		vector<vector<double> > line_mod(2);
		line_mod[0].resize(lp[0].size());
		line_mod[1].resize(lp[0].size());

		for(int ii = 0; ii < lp[0].size(); ii++){
			line_mod[0][ii] = lp[0].le1(ii);
			line_mod[1][ii] = 0.0; //MLProf.le2(ii);
			
			for (int kk = 0; kk < _Np_; ++kk) {
				
				line_mod[1][ii] += lp[kk].le2(ii);
			}
			line_mod[1][ii] /= -double(_Np_);
			//norm_fac+=line_mod[1][ii];
			
		}

		if(in_par[0][7] > 0.0){
			Convolve ff(line_mod,f2,line_mod[0].size(),f2[0].size());
			ff.xmin(line_mod[0][0]);
			ff.xmax(line_mod[0][line_mod[0].size()-1]);
			
			gsl_function F;
			F.function = &cnv_fnc;
			F.params = &ff;
			GSL_Integration itr_conv(F);
			itr_conv.xmin(ff.xmin());
			itr_conv.xmax(ff.xmax());			
			for (int i = 0; i < line_mod[0].size(); i++){
				ff.px(line_mod[0][i]);
				line_mod[1][i] = itr_conv.integrate_qag()/in_par[0][7]/sqrt(M_PI);
				itr_conv.clear();
				//cout << line_mod[0][i]
				//<< " " << line_mod[1][i] << " " 
				//<< itr_conv.integrate_qag()/2.0 << endl;
			}

		}	
				
		cout << "#Escrevendo saida em " << saida.str() << " ...\n";
		
		for (int i = 2; i < lp[0].size()-2; i++){
			int iter = (i != (lp[0].size()-1) ? 1 : -1);
			//if(lp.le2(i) > 0 || lp.le2(i+iter) > 0)
			double err = (line_mod[1][i]+cl2.flux(fase))*in_par[0][4];
			oo << line_mod[0][i]/1.0e5<< "  " << line_mod[1][i]+err*(2.0*gsl_rng_uniform(r)-1.0) //+((line_mod[1][i]+cl2.flux(fase))*in_par[0][4]*(2.*gsl_rng_uniform(r) - 1. ))
			<< " " << err << endl;
			// oo << 5000.0*(1+line_mod[0][i]/299792458.0e2) << "  " << line_mod[1][i]+((line_mod[1][i])*in_par[0][4]*(2.*gsl_rng_uniform(r) - 1. ))
			//	     << " " << (line_mod[1][i])*in_par[0][4] << endl;
			/*oo << line_mod[0][i]
			 << "  " << line_mod[1][i]+((line_mod[1][i])*in_par[0][4]*(2.*gsl_rng_uniform(r) - 1. ))
			 << " " << (line_mod[1][i])*in_par[0][4] << endl;*/
		}
		//fecha saida
		oo.close();
    }
    f_fase.close();
    }
    break;
  case 3:
    {
    DataFileReader f_star(argv[3]);
		
    surface s_star;
    while(!++f_star){
      if(f_star.numberOfColumns() != IID_NC) {
	ostringstream erro;
	erro << "Numero de colunas errado\n\tArquivo: " 
	     << argv[3] 
	     << "\n\tLido     : " << f_star.numberOfColumns()
	     << "\n\tEsperado : " << IID_NC << endl;
	aborta(erro.str());
      }
      double dado;
		int j = 1;
      s_star.theta.push_back(f_star.readColumn<double>(j));
      ++j;
      s_star.phi.push_back(f_star.readColumn<double>(j));
      ++j;
      s_star.px.push_back(f_star.readColumn<double>(j));
      ++j;
      s_star.py.push_back(f_star.readColumn<double>(j));
      ++j;
      s_star.pz.push_back(f_star.readColumn<double>(j));
      ++j;
      s_star.nx.push_back(f_star.readColumn<double>(j));
      ++j;
      s_star.ny.push_back(f_star.readColumn<double>(j));
      ++j;
      s_star.nz.push_back(f_star.readColumn<double>(j));
      ++j;
      s_star.I.push_back(f_star.readColumn<double>(j));
	  ++j;
      s_star.vx.push_back(f_star.readColumn<double>(j));
      ++j;
      s_star.vy.push_back(f_star.readColumn<double>(j));
      ++j;
      s_star.vz.push_back(f_star.readColumn<double>(j));
      //++f_star;
    }

    f_star.close();
    cout << "# Building star up with " << s_star.theta.size() << " points...\n";
    Roche cl2(s_star, in_par);
		cl2.vorb(false);
		cl2.refaz_super();
    cout << "# Done...";

    for(double fase = 0; fase <= 1; fase += 0.1){
      //prepara e abre arquivo de saida
      ostringstream saida;
      saida << argv[2] << '_' << fase;
      
      ofstream oo(saida.str().c_str());
      
      if(!oo) {
	ostringstream erro;
	erro << "Arquivo de saida nao pode ser aberto\n\tArquivo: " 
	     << argv[2];
	aborta(erro.str());
      }
      
      //escreve saida no arquivo
      
      show_par(oo, in_par);
      
      Histograma lp = cl2.rad_prof(fase, lin_form); //lin_prof;

      cout << "#Escrevendo saida em " << saida.str() << " ...\n";

      for (int i = 0; i < lp.size(); i++){
	int iter0 = ((i != (lp.size()-2) || i != (lp.size()-2)) ? 2 : 0);
	int iter1 = ((i == 0 || i == 1) ? 0 : -2);
	if(lp.le2(i) > 0 || lp.le2(i+iter0) > 0 || lp.le2(i+iter1) > 0)
	  oo << lp.le1(i)/1.0e5
	     << "  " << lp.le2(i)+in_par[1][6]+((lp.le2(i)+in_par[1][6])*in_par[0][4]*
						(2.*gsl_rng_uniform(r) - 1. ))
	     << " " << lp.le2(i)*in_par[0][4] << endl;
      }
      //fecha saida
      oo.close();
    }
    }
    break;
  case 4:
    {
    DataFileReader f_star(argv[4]);
    surface s_star;
    while(!++f_star){
      if(f_star.numberOfColumns() != IID_NC) {
	ostringstream erro;
	erro << "Numero de colunas errado em " 
	     << argv[4] 
	     << "\n\tLido     : " << f_star.numberOfColumns()
	     << "\n\tEsperado : " << IID_NC << endl;
	aborta(erro.str());
      }
      double dado;
		int j = 1;
		s_star.theta.push_back(f_star.readColumn<double>(j));
		++j;
		s_star.phi.push_back(f_star.readColumn<double>(j));
		++j;
		s_star.px.push_back(f_star.readColumn<double>(j));
		++j;
		s_star.py.push_back(f_star.readColumn<double>(j));
		++j;
		s_star.pz.push_back(f_star.readColumn<double>(j));
		++j;
		s_star.nx.push_back(f_star.readColumn<double>(j));
		++j;
		s_star.ny.push_back(f_star.readColumn<double>(j));
		++j;
		s_star.nz.push_back(f_star.readColumn<double>(j));
		++j;
		s_star.I.push_back(f_star.readColumn<double>(j));
		++j;
		s_star.vx.push_back(f_star.readColumn<double>(j));
		++j;
		s_star.vy.push_back(f_star.readColumn<double>(j));
		++j;
		s_star.vz.push_back(f_star.readColumn<double>(j));
		//   ++f_star;
    }
    f_star.close();
    Roche cl2(s_star, in_par);
	cl2.vorb(false);
	cl2.refaz_super();
	cl2.CalcPixSize();
	cout << "# - Npix = " << cl2.npix() << endl
	<< "# - Apix = " << cl2.pix_size() << endl;
		
    DataFileReader f_fase(argv[3]);
    while(!++f_fase){
      //leio fase do arquivo
      double fase;
      //f_fase >> fase;
		fase = f_fase.readColumn<double>(1);
      //preparo saida
      ostringstream saida;
      saida << argv[2] << '_' << fase;
      
      ofstream oo(saida.str().c_str());
      
      if(!oo) {
	ostringstream erro;
	erro << "Arquivo de saida nao pode ser aberto\n\tArquivo: " 
	     << argv[2];
	aborta(erro.str());
      }

      //escreve saida no arquivo

      show_par(oo, in_par);

      Histograma lp = cl2.rad_prof(fase, lin_form); //lin_prof;
      vector<vector<double> > line_mod(2);
      line_mod[0].resize(lp.size());
      line_mod[1].resize(lp.size());
      //cout << "#- Aqui 0\n";
      for(int i = 0; i < lp.size(); i++){
        line_mod[0][i] = lp.le1(i);
        line_mod[1][i] = lp.le2(i); //+flux2+in_par[1][6];
      }
      //cout << "#- Aqui 1\n";
      //tmp_vec=vector<double> (lp.col2());
      //line_mod[1].push_back(tmp_vec);
      //se tiver blurring aplica...
      if(in_par[0][7] > 0.0){
        Convolve ff(line_mod,f2,line_mod[0].size(),f2[0].size());
        ff.xmin(line_mod[0][0]);
        ff.xmax(line_mod[0][line_mod[0].size()-1]);
        gsl_function F;
        F.function = &cnv_fnc;
        F.params = &ff;
        GSL_Integration itr_conv(F);
		  
	    for (int i = 0; i < line_mod[0].size(); i++){
          itr_conv.xmin(ff.xmin());
          itr_conv.xmax(ff.xmax());
          ff.px(line_mod[0][i]);
          line_mod[1][i] = itr_conv.integrate_qag()/in_par[0][7]/sqrt(M_PI);
          itr_conv.clear();
          //cout << line_mod[0][i]
          //<< " " << line_mod[1][i] << " " 
          //<< itr_conv.integrate_qag()/2.0 << endl;
        }
	
      }
      cout << "#Escrevendo saida em " << saida.str() << " ...\n";
      for (int i = 2; i < lp.size()-2; i++){
	int iter = (i != (lp.size()-1) ? 1 : -1);
	//if(lp.le2(i) > 0 || lp.le2(i+iter) > 0) 
//	oo << line_mod[0][i]/1.0e5
//  oo << 5000.0*(1+line_mod[0][i]/299792458.0e2)
		  oo << line_mod[0][i]<< "  " << line_mod[1][i] //+((line_mod[1][i]+cl2.flux(fase))*in_par[0][4]*(2.*gsl_rng_uniform(r) - 1. ))
		  << " " << (line_mod[1][i]+cl2.flux(fase))*in_par[0][4] << endl;

	//oo << line_mod[0][i]
	 //    << "  " << line_mod[1][i]+((line_mod[1][i])*in_par[0][4]*(2.*gsl_rng_uniform(r) - 1. ))
	   //  << " " << (line_mod[1][i])*in_par[0][4] << endl;
      }
      //fecha saida
      oo.close();
    }

    }
    break;
  }
  return 0;
}

int definicoes_obscurecimento(const int law)
{
  vector<vector<double> > par;

  Light_Curve teste(par);
  
  return teste.npar_obs_borda(law);
}

void show_par(ostream& out, vector<vector<double> > &par)
{
  vector<vector<string> > par_def (4);

  par_def[0].push_back("i = ");
  par_def[0].push_back("M1 = ");
  par_def[0].push_back("M2 = ");
  par_def[0].push_back("P = ");
  par_def[0].push_back("S/N = ");
  par_def[0].push_back("Vmax = ");
  par_def[0].push_back("lres = ");
  par_def[0].push_back("iblurr = ");
  par_def[0].push_back("texp/Porb = ");
	
  par_def[1].push_back("f2 = ");
  par_def[1].push_back("grav_coef = ");
  par_def[1].push_back("limb_law = ");
  par_def[1].push_back("nspot = ");
  par_def[1].push_back("T_pole = ");
  par_def[1].push_back("ffac = ");
  par_def[1].push_back("fd = ");

  for(int i = 0; i < par[2].size(); i++){
    ostringstream str0;
    str0 << "x" << i << " = ";
    par_def[2].push_back(str0.str());
  }

  par_def[3].push_back("s_theta = ");
  par_def[3].push_back("s_phi  = ");
  par_def[3].push_back("dst = ");
  par_def[3].push_back("dsp = ");
  par_def[3].push_back("Is = ");

  out << "# Tiago Ribeiro - LinuxC++ - 2006\n#";

  for(int i = 0; i < par.size(); i++){
    for(int j = 0; j < par[i].size(); j++){
      if( i < 4 ) out << par_def[i][j];
      else out << par_def[3][j];
      
      if( i == 0 && j == 0) out << par[i][j]*180/M_PI << ' ';
      else out << par[i][j] << ' ';
    }
    out << "\n#";
  }
  out << "V(Km/s)  Flux\n";
}

int find_io(int argc, char *argv[])
{
  //Mostra Help e sai se argc < 3
  if(argc < 3) help_and_exit();
  if(argc == 3) return 1; //Modo 1: arq_in arq_out
  if(argc == 4){
    DataFileReader tmp00(argv[3]);
    ++tmp00;
    //Se numero de colunas numa linha = 1 entao argv[3] = fases.dat se for diferente considero 
    // que eh uma imagem... se nao for o programa detecta o problema na hora apropriada...
    if(tmp00.numberOfColumns() == 1) return 2; //Modo 2: arq_in arq_out fase.dat
    else return 3; //Modo 3: arq_in arq_out image.img
  }
  if(argc == 5) return 4;
  else help_and_exit();
  

}

void inline help_and_exit(void)
{
  cout << " Modo de utilizacao do programa:\n"
       << "line_prof arq_in arq_out [fase.dat] [image.img]\n"
       << "Programa para calcular o perfil de alargamento doppler das\n"
       << "linhas de absorcao/emissao de estrelas peenchendo o lobo de"
       << "Roche.\n\n"
       << "Descricao dos parametros de entrada:\n\n"
       << "arq_in o arquivo com os parametros de entrada e arq_out o\n" 
       << "arquivo para serem escritos os resultados. Opcionalmente pode\n"
       << "se incluir um arquivo com as fases desejadas para serem\n"
       << "calculados os perfis de linha ([fase.dat]). Os resultados sao\n"
       << "escritos em arquivos no formato arq_out.fase, onde arq_out o nome"
       << "\ndado na linha de entrada e fase a fase orbital. Um arquivo"
       << "arq_out\n tambem eh criado com informacoes sobre o programa.\n"
       << "Alternativamente pode-se incluir um arquivo .img (formato de imagem \n"
       << "de saida do programa roch_rom.cpp) que sera utilizado como distribuicao\n"
       << "de brilho da secundaria.\n"
       << "\nFormato do arquivo de entrada:\n\n"
       << "1a. linha de dados (sem contar linha de comentario):\n"
       << "  i    - inclinacao (graus).\n"
       << "  M1   - Massa da Primaria.\n"
       << "  M2   - Massa da Secundaria. (q = M2/M1)\n"
       << "  P    - Periodo Orbital (hrs).\n"
       << "  S/N  - Relacao Sinal-Ruido na linha.\n"
       << "  vmax - \"largura\" do histograma (km/s).\n"
       << "  lv   - \"resolucao\" do histograma (km/s).\n"
       << "  iblurr   - instrumental blurring (km/s).\n"
	   << "  texp/Porb - razao entre tempo de exposicao/periodo orbital.\n"
	   << "2a. linha de dados:\n"
       << "  f2        - fluxo da secundaria em fase 0.\n"
       << "  grav_coef - coeficiente do obscurecimento gravitacional.\n"
       << "  limb_law  - lei de obscurecimento de borda.\n"
       << "  nspot     - Numero de spots.\n"
       << "  Tpole     - Temperatura do polo da estrela.\n"
       << "  ffac      - \n"
       << "  fd        - fluxo do disco (constante somada a curva de luz).\n"
       << "3a. linha de dados:\n"
       << "  (coeficientes de obscurecimento de borda)\n"
       << "4a. a (4+nspot) linhas de dados:\n"
       << "  theta - posicao theta do spot\n"
       << "  phi   - posicao phi do spot\n"
       << "  dt    - largura a meia altura em theta\n"
       << "  dp    - largura a meia altura em phi\n"
       << "  Is    - Intensidade em relacao a I_L1\n\n"
       << "Tiago Ribeiro - LinuxC++ - 2006\n";
  exit(1);
}
